library(sna)

##load data
eventXOrg2008<-as.matrix(read.csv("Event Org Matrix 2008.csv",row.names=1,header=T,stringsAsFactors=F))
eventXOrg2009<-as.matrix(read.csv("Event Org Matrix 2009.csv",header=T,stringsAsFactors=F))
eventXOrg2009<-eventXOrg2009[,2:120]

eventAttr2008<-read.csv("Event Attributes 2008.csv",header=T,stringsAsFactors=F)
eventAttr2009<-read.csv("Event Attributes 2009.csv",header=T,stringsAsFactors=F)
eventAttr2009<-eventAttr2009[1:159,]

orgAttr<-read.csv("orgAttr.csv",stringsAsFactors=F)

##calculate one-mode projection matrices
orgXorg2008<-t(eventXOrg2008)%*%eventXOrg2008
diag(orgXorg2008)<-0
orgXorg2009<-t(eventXOrg2009)%*%eventXOrg2009
diag(orgXorg2009)<-0

##make shared-country adjacency matrix
countryMat<-matrix(0,119,119)
for(i in 1:119)
  for(j in 1:119)
    if(orgAttr$COUNTRY[i]==orgAttr$COUNTRY[j])
      countryMat[i,j]=1
diag(countryMat)<-0

##make shared-umbrella organization adjacecy matrix
umbrella<-matrix(0,119,119)
umbrella[which(orgAttr$FOE==1),which(orgAttr$FOE==1)]<-1
umbrella[which(orgAttr$GP==1),which(orgAttr$GP==1)]<-1
umbrella[which(orgAttr$WWF==1),which(orgAttr$WWF==1)]<-1
diag(umbrella)<-0

##calculate raw structural equivalence
orgXorg2008_SEDist<-sedist(orgXorg2008,method="hamming")

##makeSEL and makeSEW convert the raw structural equivalence distance matrix to the similarity matrix used by Leenders
makeSEL<-function(mat,v){
  res<-matrix(0,nrow(mat),ncol(mat))
  for(j in 1:nrow(mat)){
    s<-0
    for(i in 1:nrow(mat)){
      s<-s+(max(mat[,j])-mat[i,j])^v
    }
    for(i in 1:nrow(mat)){
      res[i,j]<-((max(mat[,j])-mat[i,j])^v)/s
    }
  }
  res
}

makeSEW<-function(l){
  res<-matrix(0,nrow(l),ncol(l))
  for(i in 1:nrow(l)){
    for(j in 1:ncol(l)){
      res[i,j]<-l[j,i]/(1-l[i,i])
    }
  }
  res
}

orgXorg2008_SEDist_w6<-makeSEW(makeSEL(orgXorg2008_SEDist,1))
diag(orgXorg2008_SEDist_w6)<-0

##normalize all matrices
orgXorg2008N<-orgXorg2008/sum(orgXorg2008)
orgXorg2008_SEDist_w6N<-orgXorg2008_SEDist_w6/sum(orgXorg2008_SEDist_w6)
countryMatN<-countryMat/sum(countryMat)
umbrellaN<-umbrella/sum(umbrella)

###make EUPER09 and UNPER09
eventNum2009<-colSums(eventXOrg2009)

targetEU09<-eventAttr2009$Target
targetEU09[which(targetEU09!=1)]<-0
targetEU09Num<-vector(length=119)
for(i in 1:119)
  targetEU09Num[i]<-targetEU09%*%eventXOrg2009[,i]

targetEUPer09<-vector(length=119)
for(i in 1:119)
  targetEUPer09[i]<-targetEU09Num[i]/eventNum2009[i]
targetEUPer09[which(is.na(targetEUPer09))]<-0

targetUN09<-eventAttr2009$Target
targetUN09[which(targetUN09!=2)]<-0
targetUN09[which(targetUN09==2)]<-1
targetUN09Num<-vector(length=119)
for(i in 1:119)
  targetUN09Num[i]<-targetUN09%*%eventXOrg2009[,i]

targetUNPer09<-vector(length=119)
for(i in 1:119)
  targetUNPer09[i]<-targetUN09Num[i]/eventNum2009[i]
targetUNPer09[which(is.na(targetUNPer09))]<-0

##get percentage of previous events that were contentious
preContPer<-vector(length=119)
preContNum<-vector(length=119)
eventNum08<-colSums(eventXOrg2008)

for(i in 1:119)
  preContNum[i]<-eventAttr2008$Contentious%*%eventXOrg2008[,i]

for(i in 1:119)
  preContPer[i]<-preContNum[i]/eventNum08[i]
preContPer[which(is.na(preContPer))]<-0


contPer<-vector(length=119)
contNum<-vector(length=119)

for(i in 1:119)
  contNum[i]<-eventAttr2009$Contentious%*%eventXOrg2009[,i]

for(i in 1:119)
  contPer[i]<-contNum[i]/eventNum2009[i]
contPer[which(is.na(contPer))]<-0


##make covariate matrices
covar<-cbind(as.numeric(orgAttr$INDMEM),as.numeric(orgAttr$STAFF),as.numeric(orgAttr$SERVICE),as.numeric(orgAttr$AGE),preContPer,targetEUPer09,targetUNPer09)
colnames(covar)<-c("IndMem","Staff","Service","Age","PrevCont","EUPer09","UNPer09")

##impute mean values for missing data in covariates
covar_noNA<-covar
for(i in 1:ncol(covar))
  covar_noNA[which(is.na(covar_noNA[,i])),i]<-mean(covar_noNA[,i],na.rm=T)

covar_int<-cbind(1,covar_noNA)

##########################
##Calculate posterior of observed data with Jeffreys prior of beta(.5,.5)
contPer_rate<-(contNum+0.5)/(colSums(eventXOrg2009)+1)
##Convert to log-odds
contPer_logit<-log(contPer_rate/(1-contPer_rate))


##run models
mod1_logit<-lm(contPer_logit~covar_int[,c(2:5,7)])
set.seed(117)
mod2_logit<-lnam(contPer_logit,covar_int[,c(1:5,7)],W1=orgXorg2008N)
set.seed(88)
mod3_logit<-lnam(contPer_logit,covar_int[,c(1:5,7)],W1=list(orgXorg2008N,orgXorg2008_SEDist_w6N))
set.seed(795)
mod4_logit<-lnam(contPer_logit,covar_int[,c(1:5,7)],W1=list(orgXorg2008N,countryMatN))
set.seed(333)
mod5_logit<-lnam(contPer_logit,covar_int[,c(1:5,7)],W1=list(orgXorg2008N,umbrellaN))
set.seed(4)
mod6_logit<-lnam(contPer_logit,covar_int[,c(1:7)],W1=list(orgXorg2008N,umbrellaN))

##presented in appendix
set.seed(387)
mod7_logit<-lnam(contPer_logit,covar_int[,c(1:6,8)],W1=list(orgXorg2008N,umbrellaN))


##### predicted probabilities
mod1_pred<-predict(mod1_logit)
mod5_pred<-vector(length=119)
for(i in 1:119){
  mod5_pred[i]<-covar_int[i,c(1:5,7)]%*%mod5_logit$beta[1:6]+mod5_logit$rho1[1]*(contPer_logit%*%orgXorg2008N[,i])+mod5_logit$rho1[2]*(contPer_logit%*%umbrellaN[,i])
}

